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We investigate the heat conductivity k of the Heisenberg spin-1/2 ladder at finite temperature 
covering the entire range of inter-chain coupling J± , by using several numerical methods and per¬ 
turbation theory within the framework of linear response. We unveil that a perturbative prediction 
k oc JJ 2 , based on simple golden-rule arguments and valid in the strict limit J± —> 0, applies to 
a remarkably wide range of J±, qualitatively and quantitatively. In the large Jx-limit, we show 
power-law scaling of opposite nature, namely, (t oc j/. Moreover, we demonstrate the weak and 
strong coupling regimes to be connected by a broad minimum, slightly below the isotropic point at 
Jl = J||. Reducing temperature T, starting from T = oo, this minimum scales as k oc T~ 2 down 
to T on the order of the exchange coupling constant. These results provide for a comprehensive 
picture of k(J±,T) of spin ladders. 

PACS numbers: 05.60.Gg, 71.27.+a, 75.10.Jm 


Introduction. Thermodynamic properties of quantum 
many-body systems are well understood, particularly in 
the vicinity of integrable points [l]. In contrast, the vast 
majority of dynamical questions in these systems remain 
a challenge to theoretical and experimental physics as 
well, in the entire range from weak to strong breaking 
of integrability. These questions consist of several timely 
and important issues such as eigenstate thermalization 
M in cold atomic gases and, as studied in this Letter, 
quantum transport and relaxation in condensed-matter 
materials. In this context, a fundamental system is the 
one-dimensional (ID) spin-1/2 Heisenberg model. It is 
relevant to the physics of quasi-ID quantum magnets flj], 
cold atoms in optical lattices Js], nanostructures 0, and 
to physical situations in a much broader context [7;, 0 ]. 

As typical for integrable systems, the energy current 
in the spin-1 /2 Heisenberg chain is a strictly conserved 
quantity 0 [To[. This implies purely ballistic flow of heat 
at any temperature and provides the theoretical basis 
for explaining the colossal magnetic heat conduction ob¬ 
served experimentally in quasi-ID cuprates HM1 . In 
contrast to heat flow, spin dynamics, including the exis¬ 
tence of ballistic }15l-l27| and diffusive transport channels 
|, is theoretically resolved only partially, and also 


under ongoing experimental scrutiny [34 

Because of strict energy-current conservation in this 
model, the heat conductivity n is highly susceptible to 
breaking of integrability by, e.g., spin-ph onon coupling 
[Ud , dimerization or disorder , and interac¬ 

tions between further neighbors [4a, |4£j|. One °f the 
most important perturbations is inter-chain coupling, i.e. 
Jl, which is the key ingredient to spin-ladder compounds 
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FIG. 1. (Color online) Thermal conductivity k (per chain) 
versus J±/J\\ and /?. PT: known perturbative regime. Issues 
clarified in this Letter: extent of power-law scaling (dashed 
line) close to PT and sum-rule (SR) regimes, nonperturbative 
numerical treatment for entire Jx-range, location of minimal 
conductivity, and temperature variation. 


0,0- Since the discovery of the discontinuous tran¬ 
sition from one to two dimensions in quantum magnets 
d, spin ladders are a cornerstone of correlated elec¬ 
tron systems. They display quantum confinement ptij . 
transforming gapless spinons of simple spin chains into 
new massive triplons [49|, |50( ■ They provide insights into 
fractionalization, quantum phase transitions 51], Bose- 
Einstein condensation [521 . and disorder-induced mag¬ 
netism !53|. They are paradigmatic to high-Tc super¬ 
conductors, undoped [54| and doped j,55]]- They serve as 
models in other fields, e.g., cold atomic gases [56|, quan¬ 
tum information theory [57j , and carbon nanotubes [Hit . 

Early on, perturbation theory (PT) to lowest order, 
i.e., a simple golden-rule argument |59|,[60j], has suggested 
dissipative heat flow with a scaling k oc JJ 2 , as illustrated 
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on the l.h.s. of Fig. [T] However, the relevance of such scal¬ 
ing is unclear off the strict limit J± —> 0, as is the radius 
of convergence of the PT. Understanding k over a wider 
J± range has been hampered by the lack of sufficiently 
accurate nonperturbative methods. In particular, state- 
of-the-art numerical methods have been restricted to the 
regime J± = 0(1), where finite-size effects are moder¬ 
ate and spectral structures are broad, i.e., time scales 
are short [6jJ, |62|. Thus, heat transport in the transition 
from weakly coupled chains to strongly coupled ladders 
is understood only in few and narrow regions. 

In this Letter, we lift these restrictions and study the 
heat conductivity n over the entire range of the inter¬ 
chain coupling J±. Using several methods within lin¬ 
ear response, we (a) quantitatively connect to PT in the 
small- Jj_ limit and (b) unveil its validity for a remark¬ 
ably wide range of J±. In addition to the PT, scaling 
as k oc JJ 2 , we (c) demonstrate a qualitatively different 
power-law scaling k oc Jj_ in the large Jj_-limit. Conse¬ 
quently, we (d) find a broad minimum of n in the region 
J_l < 1. Reducing temperature T, starting from T = oo, 
this minimum (e) scales as k oc T~ 2 down to T on the 
order of the exchange coupling. Thus, we provide a com¬ 
prehensive picture of n(J±,T ), beyond the known results 
sketched as part of Fig. Q] 

Model. We study a Heisenberg spin-1/2 ladder of 
length TV/2 with periodic boundary conditions. The 
Hamiltonian H = TLy + H± consists of a leg part 7/ 
and a rung part H±, 

z n/2 N/2 

H w = J n Si ’ k ■ , h ± = j±J2 s o - Si ’ 2 ’ w 

k— 1 i —1 i =1 

where S^ are spin-1/2 operators at site ( i,k ), J || > 0 
is the antiferromagnetic leg coupling, and J± > 0 is the 
rung interaction. 2 = 2 is the number of legs. For Jj_ = 0, 
the ladder splits into integrable chains, with a gapless 
ground state and spinon excitations. For J|| = 0, it sim¬ 
plifies to uncoupled dimers, with a gapped ground state 
and triplon excitations. For Jj_, Ju ^ 0, the ladder is non- 
integrable. Generally, the model in Eq. (HD preserves the 
total magnetization S z and is translationally invariant. 
We focus on the representative sector S z = 0 [?tI |. 

The energy current has the well-known form j = y +j± 

fill, 


2 N/2 

j \| = J | 'y [ y [ • (Sx Si+qfc), (2) 

k= 1 2=1 

J J z N/2 

J± = EE( S -U - Si+i,*) • (Sq fc X S,3- fc ) • 

fc=1 i=l 

j and H commute only at the integrable point J± = 
0. We investigate the autocorrelation function at inverse 


temperatures /? = 1/T, 


(3) 


TV Tr{e-P H } 


where the time argument of j ( t ) refers to the Heisenberg 
picture, j = j( 0), and C(0) = 3(Jj| + J|J 2 /2)/32 for 
pj\\ —> 0. 

From C(t ), we hrst determine the Fourier transform 
C(u>) and then the conductivity via the low-frequency 
limit k/z = /3 2 C(oj 0). Additionally, we can extract 
the conductivity directly by k/z = /3 2 f * 1 d tC(t). Here, 
the cut-off time t\ has to be chosen much larger than the 
relaxation time r, where C'(r)/C(0) = 1/e [631 ]. 

Methods. We calculate C by complementary numerical 
methods, with a particular focus on dynamical quantum 
typicality (DQT) [25, 2^,[64| (see also Refs. l65l - l73h . DQT 
relies on the time-domain relation 


C(t) = Re 


(<l>/i(t)\jWp(t)) 
IV(^(0)|^(0)) + ’ 

>, I Vp{t)) = e~ lHt je-P H / 2 


(4) 


| <f>f}(t)) = e - lHt ~P H / 2 \ 
where | ip) is a single pure state drawn at random and 
e scales inversely with the partition function, i.e., e is 
exponentially small in the number of thermally occu¬ 
pied eigenstates 2^, [26, 64j. The great advantage of Eq. 
© is that it can be calculated without any diagonaliza- 
tion by the use of forward-iterator algorithms. We use a 
fourth-order Runge-Kutta iterator with a discrete time 
step St J|| = 0.01 <C 1. Together with sparse-matrix rep¬ 
resentations of operators, we can reach systems sizes as 
large as N = 32. For more details on the method and its 
accuracy, see Refs. 0 andfTTl 

Additionally, we confirm our DQT results with numeri¬ 
cal methods based on Lanczos diagonalization in the fre¬ 
quency domain [76[ . with the frequency resolution Sui 
crucially depending on the number of Lanczos steps M, 
Suj oc 1/M. At low T, we choose the finite-T Lanczos 
method (FTLM) with M ~ 200 771. At high T, we also 
use the microcanonical Lanczos method (MCLM) with 
M ~ 2000, significantly improving Slo. 

Results. We begin with J±/J\\ > 1 and /3Jy —► 0. In 
Fig. [2] (a) we summarize our DQT results on C(t) for 
different J±/J\\ = 1, 1-5, 2. Several comments are in 
order. First, the initial value C(0) agrees with the high- 
T sum rule and therefore increases with J±. Second, all 
C(t) depicted decay to zero on a time scale 5r ~ 10/ Jy. 
Third, the C(t) curves do not change when the number of 
sites is increased from TV = 22 to 32. Thus, we observe 
very little finite-size effects, i.e., we can safely consider 
our results as results on C(t) for TV —>• oo. Note that for 
TV > 30 we consider a single translation subspace k since, 
for these TV, C(t) is k independent at /3 —> 0 mm. 

Next, we discuss the spectrum C(co). To this end, we 
show in Fig. [2] (b) for J±/J || = 1, 2 the Fourier transform 
of our DQT data for times t < lOr ~ 20/J||. These 















3 


K 5 



4 6 

Time tJ\\ 


10 


1.2 


0.8 


3 

O 0.4 


N = 22 (all k) 0 3 \ 

0 - N = 28 Ref. 49 

" ■ N = 30 (fc = 0) 

V* N = 32 (k = 0) a 


Jl 

,N _ --10-.. 

/Jy = 1,2 0 0-3 0.6 

... (b) a/-domain • 

, 0 «o ♦ o» O «C> 


0.5 1 1.5 

Frequency u /./ 


- N = 22 (all k) 

- N = 30 (k = 0) 

- Jl — 0 Sum rule 



O 


J±/J\\ = 0.75,0.5,0.25,0.15 
(a) 1-domain 


50 100 

Time tJn 


150 


2.4 

1.8 

1.2 

0.6 

0 


N = 22 (all k) 2 
■-*- N = 28 (all k) 1 ' 
♦ N = 32 (k = 0) 1 



MCLM 


**Xxxwy» 

' 0 0.1 0.2 
(b) tj-domain 
JaJJ\\ = 0.25 


0.1 0.2 0.3 

Frequency ui/ Ju 


0.4 


FIG. 2. (Color online) The (a) t and (b) u dependence of the 
autocorrelation C for strong J±/J\\ > 1, /3J\\ -4 0, and N < 
32, as obtained from DQT. Spectra in (b) are obtained by 
Fourier transforming finite-f data t < lOr ~ 20/Jy (symbols); 
Inset: Low-ui limit for J±/J\\ = 1, the largest N = 32, and 
t < 5r (crosses), lOr (other symbols), 50 t (curves); Main 
panel: Spectra from Lanczos methods (N = 28 MCLM of Ref. 
[6l] , N = 22 FTLM) are shown (curves). Note that method- 
related errors are negligibly small 0- 


FIG. 3. (Color online) (a) The t dependence of C for vari¬ 
ous small J±/J || = 0.15, ..., 0.75, obtained from DQT for 
PJ\\ —¥ 0 and N < 30. (b) Spectrum for J±/J\\ = 0.25, ob¬ 
tained by Fourier transforming fmite-t data t < 5r ~ 80/Jy 
(symbols); Inset: Low-cj limit for the largest IV = 32 and 
t < 5t (diamonds), lOr (crosses); Main panel: Spectrum 
from N = 22 and 28 MCLM and a Lorentzian fit are shown 
(curves). 


times correspond to a frequency resolution Su> ~ 0.15J||. 
For this resolution, the Fourier transform is a smooth 
function of u> and displays a well-behaved limit for io —> 
0, i.e., C(lo — > 0) = C(oj = 0). Moreover, this limit 
and C(ui) in general do not depend on system size for 
N > 22. The inset of Fig. [2] (b) clarifies the impact of the 
lo resolution by displaying additional Fourier transforms 
of DQT data, evaluated for shorter (t < 5r) and longer 
(t < 50t) times at J±/J\\ = 1 and for the largest N = 32. 
Clearly, the low-cu limit is independent of the to resolution 
resulting from the specific choice of t. This robustness, 
together with the N independence, allows us to reliably 
extract a quantitative value for the dc conductivity at 
Jl/J|| = 1, K/z/3 2 Jf = 0.29. 


To additionally demonstrate the validity of our DQT 
approach, we compare to our FTLM results and to exist¬ 
ing MCLM spectra from the literature 61] in Fig. [2] (b). 
Obviously, the agreement is very good. 


Now, we turn to small J±/J\\ < 1. In Fig. [3] (a) we 
depict our DQT results on C[t) for various J±/J\\ = 
0.15,..., 0.75. The initial value C(0) approaches the 
J_l = 0 sum rule when J± is reduced. Furthermore, the 
decay is slower for smaller J± and finite-size effects are 
naturally stronger in the vicinity of the integrable point 
Jl = 0 . For the smallest J±/J\\ = 0.15 depicted, these 
finite-size effects are still moderate when comparing C[t) 
for N = 22, 30. In Fig. [3] (b) we show the Fourier trans¬ 


form of C(t < 5t ~ 80J||) for J±/J\\ = 0.25. For the 
largest N = 32, this Fourier transform is well described 
by a Lorentzian line shape and, again, the low-w limit 
does not depend on t. Since C(oj) has a narrow spec¬ 
trum, MCLM with a high ui resolution (M = 2000) is a 
better choice for comparison than FTLM (M = 200) a 
and agrees well with DQT. Note that resolving narrow 
spectral features by DQT is a new concept of our Letter, 
which can be applied in a much broader context. 

Next, we discuss the scaling of the conductivity n over 
the entire range of J±. In Fig.[|](a) we summarize k(J_l) 
as inferred from DQT data for C(t < 5r). Here, we ob¬ 
serve a broad minimum of k(Jj_), centered between two 
regimes with power-law scaling at large and small J±. 
The scaling oc J\ in the large J± limit is a direct conse¬ 
quence of the static sum rule C(0) oc J\, noted following 
Eq. ©. The scaling oc JJ 2 for small J_l, however, is not 
simply related to (7(0) since (7(0) « const, for such J±. 
Particularly, we find this scaling to hold over a remark¬ 
ably wide range of 0.07 < J±/J\\ < 0.35. This finding is 
a central result of this Letter. Below J±/J|| < 0.07, com¬ 
putational efforts for 5r data are very high and finite-size 
effects are too large, even for N accessible to DQT. 


To gain further insight into the scaling at small J_l, we 
calculate the scattering rate 7 = 1 /r to lowest order in 
J_l, i.e., J 2 , following the PTs in Refs. [5§, 6 (| 74, and 
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FIG. 4. (Color online) (a) Scaling of the conductivity k with 
J_l, obtained from DQT and finite-f data t < 5r for /3J m —» 0 
and JV < 32 (closed symbols). Results for the simplified oper¬ 
ator j' — j || are also depicted at J±/J\\ ~ 1 (open symbols). 
Additionally, power laws 0.097(Jx/./||)~ 2 and 0.21(Jx/J||) 2 
are shown (lines), (b) PT for the scattering rate 7 , carried 
out using DQT. The PT of Ref. l59i is also depicted (bullet). 


751 . This rate reads (/3J\\ —t 0) 


7 = lint 

1 1 —>-00 


d t 


Tr{z[j||, i7x](I||) H±\} 

Tr{.lU- 


oc J 


_L > 


(5) 


where 1 1 | refers to the Heisenberg picture of H ||. Figure 
3] (b) shows 7 evaluated by DQT applied to Eq. ([5]) for 
large N < 30. Note that this application of DQT is a new 
concept of our Letter 77]. As shown in Fig. [2 we find 
good agreement with previous evaluation of 7 in Ref. [5fj 
based on smaller systems. Most notably, however, 7 well 
agrees with the scattering rate 7 ' as extracted directly 
from k in Fig. 2] (a) via the relation 7 ' = z/3 2 C(0 )/k. 
This agreement is another main result of our Letter. Note 
that PT holds up to J±/J\\ ~ 1 for the simplified current 
j = j ||, see Fig. [4] (a), which is the regime where the 
system behaves Markovian, i.e., has no memory. For the 
explicit calculation of the memory kernel, see Ref. tzJ 
Now we turn to /3J\\ ^ 0, focusing on J±/J\\ = 1. 
In Fig. [5] (a) we depict our DQT results for C[t) for 
/3J|I = 0.5,..., 1.5. While C(0) decreases as /3 is in¬ 
creased, the relaxation time shows the tendency to in¬ 
crease with /3. However, significant finite-size effects ap¬ 
pear as nondecaying Drude weights. Since these Drude 
weights exceed 20% of (7(0) at 3J\\ ~ 1.5, we restrict our¬ 
selves to /3J\\ < 1. For such /?, once again, FTLM agrees 
with the Fourier transform of our DQT data, which also 
shows a N independent dc limit for large N ~ 30, see 
Fig. [5] (b). Finally, in the inset of Fig. [5] (a) we show 
the T dependence of the conductivity k. Remarkably, 
in the T range accessible to our methods, we observe no 
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FIG. 5. (Color online) (a) The t dependence of C for /3J\\ = 
0.5,..., 1.5, obtained from DQT for J±/J\\ = 1 and N = 28. 
(b) Spectrum for fiJ\\ = L obtained by Fourier transforming 
finite-t data t < 5r ~ 12/Jy for N < 32. Additionally, a 
spectrum from N = 22 FTLM is depicted (curve). Inset: T 
dependence of the conductivity k, calculated by N = 32 DQT 
(closed symbols, curve) and N = 22 FTLM (open symbols). 


significant deviations from the high-T behavior n oc /3 2 . 
While these T are low from a numerical point of view, 
they are still too high for a comparison to experiments 
on yet available materials, where the exchange coupling 
constant is large. 

Conclusion. We studied the heat conductivity k of 
the Heisenberg spin-1/2 ladder at finite temperature and 
over the entire range of the rung interaction Jj_, using 
several methods within linear response. We detailed the 
power-law scalings n oc JJ 2 and k oc J| at weak and 
strong J l, respectively. We found a broad minimum of 
k in the region J± ~ 1, with a scaling of its temperature 
dependence as k oc T~ 2 down to T on the order of the 
exchange coupling. Thus, we provided a comprehensive 
picture of k(J±,T). 
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SUPPLEMENTAL MATERIAL 

PERTURBATION THEORY 
Leading-Order Scattering Rate in the Markov Limit 


Note that, in the Markov approximation, the qualitative 
scaling oc 1 / Jj_ does not depend on details of the memory 
kernel K(t) while the quantitative value of the scattering 
rate 7 clearly does. 


In this section we discuss the perturbation theory for 
the energy current j = ju + j± in detail. In the limit of 
small inter-chain coupling, J± —> 0, the leg part j || is the 
dominant contribution, i.e., 

j=3\\+0(J±). (SI) 

Because the leg part j\\ is strictly conserved for the leg 
Hamiltonian H \|, [j||,H)|] = 0, the rung Hamiltonian H± 
is the only origin of the scattering of j\\. This scattering 
can be treated perturbatively if the inter-chain coupling 
Jj_ is a sufficiently small parameter. In the time domain, 
we can formulate such a perturbation theory in terms of 
the integro-differential equation 

t 

C(t) = - Jdt' K(t - t') C(t) (S2) 

0 


Strong Perturbations 

For strong inter-chain coupling J±, the perturbation 
theory discussed above necessarily breaks down for two 
different reasons: (i) The memory kernel K(t) in Eq. (IS3I) 
does not incorporate higher-order contributions, (ii) The 
current operator j is approximated by the leg part j \\. In 
fact, 

j ~ j -l (S9) 

in the limit of large J±. In this limit, C(0) = const, turns 
into C(0) oc Jj_. This scaling with J± reflects that the 
dominant energy contribution is the bond energy in the 
rungs. 

Dynamical Quantum Typicality 


for the autocorrelation function C(t) = Re(j|| (t)j\\)/N of 
the leg part j \\, where the memory kernel K(t) occurs in 
the time convolution on the right side. To leading order 
of the perturbation J±, J^_, and in the high-temperature 
limit, f3 —> 0, this memory kernel reads fsiil . 01 


K(t) 


Tr{»[j|| 1 ) i[j \\, H±]} 

Tr {j|} 


oc J 2 


_L ! 


(S3) 


where 1 1 | indicates the Heisenberg picture with respect to 
H ||, i.e., 


Assuming that K{t) fully decays on a finite time scale tk, 
i.e., using the Markov approximation, Eq. (IS2I) simplifies 
for small J ±, where C(t) decays on a very long time scale 
r > t K - Thus, Eq. (|52]> factorizes into 

C(t) = - 7 C(t), (S5) 


where 7 is the scattering rate 
1 ft 

7 = — = lim / d t' K[t’) oc , (S6) 

T t—^OO 7 Q 

implies the 

(S7) 


cf. Eq. (5) of our Letter. Obviously, Eq. 
exponential relaxation 

C(t) = C(0) . 

Consequently, the heat conductivity becomes 


-752 = lim [ dt' C(t') 

Z/3 2 t^ooj 0 


— “ 4 - (S8) 


Verifying the Markov approximation and determining 
the quantitative value of the scattering rate 7 necessarily 
requires full knowledge about the time dependence of the 
memory kernel K (t). Even though this time dependence 
is generated by the integrable Hamiltonian H ||, cf. Eqs. 
(IS3l) and (|S4|) . an exact calculation of K[t) is unfeasible 
due to the complexity of H\\. Therefore, in praxis, K(t ) 
has to be calculated numerically. The standard approach 
is the exact diagonalization of [si, S2]. However, since 
H || is a many-body Hamiltonian, this approach is only 
feasible for at most IV ~ 16 sites and finite-size effects 
can be large for such N. 

To overcome the limitation of exact diagonalization to 
small IV, we first need to note that the memory kernel 
K(t) for /3 —0 is just the autocorrelation function of the 
Hermitian operator 


V Tr{j 1 } ’ 


(S10) 


i.e., K(t) = Tr{j'(t||)j'}. Hence, remarkably, we can use 
the concept of dynamical quantum typicality to calculate 
K(t). Specifically, in analogy to Eq. (4) of our Letter, we 
get the relation for /3 —> 0 


K(t) 


(m\rwt)) 

(m\m) 


with the two auxiliary states 

kW)=e-^"‘/|V’), 


(Sll) 


7 


(512) 

(513) 
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FIG. SI. (Color online) Full time dependence of the memory 
kernel K(t), according to the leading-order prediction in Eq. 
(OH, for high temperatures 0 —> 0 and different lattice sites 
N = 22, 26, 30. The data depicted is numerically calculated 
using the typicality relation in Eqs. EH i, EH, EH and a 
fourth-order Runge-Kutta iterator with a discrete time step 
5tJ || = 0.01 -C 1. Apparently, K(t) does not depend on N 
and fully decays on a rather short time scale tk J\\ ~ 2 . As 
compared to exact diagonalization of, say, N = 16 sites, the 
Hilbert-space dimension accessible is larger by a factor of ca. 
16,000. 


where \if) is a single pure state drawn at random. Again, 
e scales inversely with the partition function, i.e., e is 
exponentially small in the number of thermally occupied 
eigenstates S5-|S7|. 

The typicality relation in Eqs. (IS11D . (IS12I) . and (IS13I) 
can be calculated for as many sites as N = 30, using a 
fourth-order Runge-Kutta iterator with a discrete time 
step St J|| = 0.01 -C 1 and sparse-matrix representations 
of the operators H\\ and j' [S7J. In Fig. [Sl]we depict our 
results on the time dependence of K(t ) for 0 —> 0 and 
different N. Apparently, K (t) does not depend on N and 
fully decays on a rather short time scale tk J\\ ~ 2 . Thus, 
the Markov approximation is indeed justified. Note that 
the area under the K(t) curve is the scattering rate 7 
shown in Fig. 4 (b) of our Letter. 
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FIG. S2. (Color online) Error analysis: The time-dependent 
energy-current autocorrelation function C(t) is numerically 
calculated according to the typicality relation in Eq. (4) of 
our Letter for two random initial states 7 ^ |i/j) 2, N = 30 
lattice sites, equal exchange couplings J± = Jy , and different 
temperatures (a) 0 Jy = 0, (b) /5 Jy = 0.5, (c) 0 J\\ = 1. In 
addition, results for N = 32 and a single initial state are 
shown. All results in (a)-(c) correspond to the sector S z = 0 
and k = 0. Insets: Relative deviation of the initial value C(0), 
as obtained from \ip )2 or N = 32, from the C(0) of |-0) 1. (Note 
that the x axis is meaningless.) For the 0 range depicted, the 
initial-state dependence is much below 1% and smaller than 
the, also small, finite-size effects. 


ERROR ANALYSIS 


p(e) is upper bounded by |S5HS7| 


Specific Realization of the Initial State 

Our main numerical method used is essentially based 
on the typicality relation in Eq. (4) of our Letter, where 
the random error e occurs. In this section we discuss this 
error in detail. The probability distribution p(e), i.e., the 
probability to get an error of size e, has the mean 

e = 0 . (S14) 

Hence, if averaging is performed over sufficiently many 
random initial states | ip)i, then any error vanishes. Note 
that we do not need to perform such averaging for reasons 
outlined below. The width of the probability distribution 


<r(e) < O 


( \/Rc (j(t)jj(t)j) \ 

\ N y/dtfi ) ’ 


where the effective dimension 


d eS = Tr{e-^ H ~ E ^} 


(S15) 


(S16) 


is the partition function and Eq denotes the ground-state 
energy. Hence, the maximum error a(e) decreases faster 
with system size than l/y/d c s does. In the limit of high 
temperatures, 0 —> 0 , d e s = 2 W is a huge number for 
TV ~ 30 and the maximum error consequently is tiny. In 
fact, it has been shown in Ref. |s3 that significant errors 
only occur for effective dimensions below d e s ~ 10 , 000 , 







































e.g., for rather small N. However, for large N, d e g can 
also become small for two reasons relevant to the study 
in our Letter. 

(i) To reduce computational effort for large N > 30, 
we restrict our investigation to a single but representative 
symmetry subspace, i.e., to the quantum numbers S z = 0 
and k = 0. While this restriction does not have impact on 
the exact autocorrelation function C(t ), the dimension of 
the subspace 


do,o^^J«2" 


N/2 \N/2 


(S17) 


is much smaller than the full Hilbert-space dimension. In 
the high-temperature limit, /3 —> 0, d e ff = do.o is still a 
large number for N ~ 30. 

(ii) If temperature is reduced from infinity at fixed N, 
d e ff gradually deceases and eventually becomes 1 at zero 
temperature. Therefore, d e s becomes a small number for 
sufficiently low temperatures and, as a consequence, the 
upper bound in Eq. (IS15I) does not imply a small e any 
further. Note that e can still be small since e = 0 at zero 
temperature. 

Due to (i) and (ii), it is important to verify in praxis 
that e is indeed a negligibly small error. This verification 
is most conveniently done by repeating the calculation 
of the autocorrelation function C(t) for a second random 
initial state |i />)2 ^ | ip)i- In Fig. [S2]we depict C(t), as 
obtained from and \ip) 2 , for N = 30 lattice sites, 
quantum numbers S z = 0 and k = 0, equal couplings 
J± = J||, and different temperatures (3 J|| < 1. For this 
temperature range, we extract the heat conductivity in 
our Letter. Clearly, the C(t) curves for \4>)i and \4>)2 in 
Fig.[S2]are very close to each other. This independence of 
the specific realization of the random initial state proves 
a small e for temperatures /3 J|j = 1. 

It is also instructive to quantify the size of errors. To 
this end, let us consider the relative error of the initial 
value given by 


Cr(0) 


|C(0,|V.) 2 )-C(0,M 1 )| 

c( o,|V>}i) 


(S18) 


As shown in the insets of Fig. IS21 e r (0) is much smaller 
than 1% and particularly smaller than the, also small, 
finite-size effects. 


Runge-Kutta Time Step 

The typicality relation in Eq. (4) of our Letter requires 
to propagate pure states in real and imaginary time. We 
perform the propagation by a fourth-order Runge-Kutta 
iterator with a discrete time step St J\\ = 0.01 <C 1. This 
time step is a potential source for errors if the relaxation 
time of the energy current is very long, i.e., for very small 
inter-chain couplings J±/J\\ -C 1. To ensure sufficiently 


high accuracy, we verified for all J± that the norm of the 
two pure states and \</>p(t)) propagated does not 

deviate significantly from 1 on times up to the relaxation 
time of the energy current. In Fig.[S2]we illustrate that, 
even for the demanding case J±/J\\ = 0.1, we find that 
deviations from 1 are less than 1%. Note that reducing 
the time step for the large N = 30 depicted is unfeasible 
for J±/J\\ ~ 0.1. 


Lanczos-Related Errors 


Within the Lanczos-diagonalization techniques used in 
our Letter, the origin of potential errors is twofold and 
related to [S8I |: (i) spectral resolution Sui and (ii) number 
of effective terms in the thermodynamic sum Z*. 

(i) The spectral resolution is given by 


5u> = 


A E 

~M 


(S19) 


where M is the number of Lanczos steps used. Note that 
we use M = 200 for FTLM and M = 2000 for MCLM in 
our Letter. A E = E max — E m i n is the full energy span of 
the Hamiltonian, i.e., E m i n and E max are the smallest and 
largest eigenenergy, respectively. This span depends only 
weakly on the inter-chain coupling Jj_, e.g., for N = 22 
we find AE/J\\ ~ 30 for the ladder case J±/J\\ = 1 and 



0 200 400 600 800 1000 

Time tJ« 

FIG. S3. (Color online) Error analysis: Real-time dependence 
of (a) energy-current autocorrelation function C(t) and (b) 
norm a/(4 >/3(t)|4>/3(t)) of the pure state |$^(t)) in Eq. (4) of 
our Letter for small inter-chain coupling J±/J\\ = 0.1, N = 30 
lattice sites, sector S z = 0 and k = 0, and high temperatures 
P —> 0. Although C(t) decays on a long time scale, the norm 
\J (®^(t)|4 > ^(t)) does not deviate more than 0.25% from 1 on 
this time scale. Hence, our choice of the Runge-Kutta time 
step 5t J|| = 0.01 <§; 1 is reasonable. 
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FIG. S4. (Color online) Error analysis: FTLM results for the 
frequency-dependent energy-current autocorrelation function 
C(u>) for different Lanczos steps M = 50, 100, 150, 200 and 
two temperatures (a) /3 Ju =0 and (b) P J|| = 1. All results 
depicted correspond to strong coupling J±/J\\ = 1, N = 22 
lattice sites, and the sector S z = 0 (all k). 


AE/J\\ « 20 for the chain case J±/J\\ = 0. This span, 
together with M = 200, yield Soj/J^ ss 0.15. It is evident 
from Fig. 2 of our Letter that this spectral resolution is 
sufficient for large J±/J\\ = 0(1). However, the spectrum 
C(oj) is narrow for small Jj_/J\\ = 0.25 and has a width 
of roughly 0.2 J\\. Thus, we have to turn to MCLM and 
M = 2000 for sufficiently high spectral resolution. 

In Fig. [S4] we show the M dependence of the FTLM 
result for strong coupling J±/J\\ = 1, two temperatures 
P J|| = 0, 1, and N = 22 lattice sites. Obviously, finite M 
are visible as quasi-periodic oscillations. But the overall 
structure of C(oj) and the dc limit C(oj —> 0) are already 
well converged for M = 200. 

(b) Similar to the typicality approach discussed before, 
the statistical error of the Lanczos procedure is 

(S20) 

where R is the number of random pure states used for 
sampling and Z* is the thermodynamic sum 

Z* = ^2e- 0iEn ~ Eo) , (S21) 

n 

where E 0 is the ground-state energy. Note that, for /3 = 0 
and /3 = oo, the thermodynamic sum Z* is equivalent to 
the effective dimension d e s in Eq. (1S16I) . For any finite 
/ 3 , however, Z* « d e s since energies E n located in the 
middle of the spectrum are approximately correct within 
the Lanczos technique. We find Z* ss 450 for Jj_/Jn = 1, 


N = 22, and /3 J|| = 1 and sample over R = 10 random 
pure states for the /? yf 0 cases in Fig. 5 our Letter. For 
all /3 = 0 cases, R = 1. 


In Fig. [S5] we depict FTLM results for two different 
random states |$)i ^ |$) 2 (R = 1) and an average over 
several |$)j (R = 10) for strong coupling J±/J\\ = 1, two 
temperatures /3 J|| =0, 1, and N = 22 lattice sites. It is 
evident that the dependence on R is negligibly small. 

For a more detailed description of the implementation 
of FTLM and MCLM, we refer the interested reader to 
Refs. Issl-lsTfi 


EQUIVALENCE OF ENSEMBLES 

Finally, we also demonstrate that all results presented 
in our Letter do not depend on the restriction to the 
magnetization sector S z = 0 (canonical ensemble). To 
this end, we repeat the calculation in Fig. 2 (b) of our 
Letter for ( S z ) = 0 (grand-canonical ensemble), taking 
into account all magnetization sectors. The result of this 
calculation is depicted in Fig.lSGland proves that S z = 0 
and ( S z ) = 0 yield the same physics. 

For the J\\/J± = 1 and N = 30 DQT spectrum shown 
in Fig. [SH we also depict in Fig. [57] (a) the underlying 
real-time data. Furthermore, we comp are this real-time 
data to the tDMRG data of Ref. [Sill and find excellent 
agreement at /3J|| = 0. As illustrated in Fig. [57] (b), the 
agreement between DQT and tDMRG data is also very 




0 0.5 1 1.5 2 

Frequency u/ J\\ 

FIG. S5. (Color online) Error analysis: FTLM results for the 
frequency-dependent energy-current autocorrelation function 
C(w) for two different random states |^>)i ^ |$)2 (R = 1), an 
average over several |<I>); (R = 10), and two temperatures (a) 
P Jn =0 and (b) P Ju =1. All results depicted correspond 
to strong coupling J±/J\\ = 1, N = 22 lattice sites, and the 
sector S z = 0 (all k). 
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FIG. S6. (Color online) The same as Fig. 2 (b) of our Letter 
but for the grand-canonical ensemble ( S z ) = 0. 


good for /3J|| = 1. Recall that the Fourier transforms in 
the inset of Fig. [S6l rely on much longer times than those 
times depicted in Fig. 1571 


TEMPERATURE DEPENDENCE 

To illustrate that the temperature dependence of the 
heat conductivity does not depend on the specific point 
J±/J\\ = 1 considered in our Letter, we also repeat the 
calculation in Fig. 5 for J±/J\\ = 0.5. The results of this 
calculation are shown in Fig. IS8I Most importantly, we 
again find the scaling n ex T -2 down to T on the order 
of the exchange coupling, as shown in the inset of Fig. 
EH (a). 



0 2 4 6 8 10 

Time tJii 

FIG. S7. (Color online) (a) Real-time data underlying the 
J±/J\\ = 1 and N = 30 DQT spectrum in Fig. IS6I (b) The 
same as (a) but for /3Jy =1 ^ 0. For comparison, in (a) and 
(b) the tDMRG data of Ref. ISllI is depicted. 



0 12 3 

Frequency to/ J\\ 

FIG. S8. (Color online) Real-time decay of the energy-current 
autocorrelation function C(t) for /3Jy = 0,0.5,1, obtained 
from DQT for = 0.5 and N = 30. (b) Spectrum for 

/3J || = 1, obtained by Fourier transforming data for finite 
times t < 5t ~ 30/J||. Inset: Temperature dependence of 
the conductivity k, calculated by N = 30 DQT. The overall 
figure is similar to Fig. 5 of our Letter, where J±/J\\ = 1. 
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